#include "header.h"

/* -------------------------------------------------------------------------- */
void Line(){
/* -------------------------------------------------------------------------- */

printf("-------------------------------------------------------------------\n");

}

/* -------------------------------------------------------------------------- */
void title(){
/* -------------------------------------------------------------------------- */

Line();
printf("                        3D KMM Code for External Flows            \n");
Line();
printf("\n");
fflush(NULL);

}
/* -------------------------------------------------------------------------- */
void output(struct node* p){
/* -------------------------------------------------------------------------- */
int i,point;
double fact,totemp,totemp_inf,totpres,totpres_inf,pres;
double rho,temp,u,v,w,mach,sonic,vel,entropy,entropy_inf;
double q2_inf, cp;
FILE *fpt, *fsolid, *fcp;

fpt = fopen(OutFile,"w");
fsolid = fopen("solid.dat", "w");

fprintf(fsolid, "Variables=X Y Z Pressure Mach u  v  w\n");
for(i=0;i<no_grid_points;i++){
	pres  = p[i].var[4];
	fprintf(fpt,"%20.10e %20.10e %20.10e %20.10e %20.10e %20.10e %20.10e %20.10e\n",
			p[i].x,p[i].y,p[i].z,p[i].var[0],p[i].var[1],p[i].var[2],p[i].var[3],pres);
	if(p[i].point_type==SOLID_WALL || p[i].point_type==KUTTA){
		u = p[i].var[1];
		v = p[i].var[2];
		w = p[i].var[3];
		vel = sqrt(u*u + v*v + w*w);
		mach = vel/sqrt(GAMMA*pres/p[i].var[0]);
		pres = (pres - pre_inf)/q_inf;
		fprintf(fsolid,"%15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e %15.6e\n", 
					 p[i].x, p[i].y, p[i].z, pres, mach, u, v, w);
		}
	}
fclose(fpt);
fclose(fsolid);

fpt = fopen(PointFile,"w");

fact = 1.0 + 0.5*(GAMMA-1.0)*mach_inf*mach_inf;
totemp_inf  = temp_inf*fact;
fact = pow(fact,GAMMA/(GAMMA-1.0));
totpres_inf = pre_inf*fact;
entropy_inf = pre_inf/pow(rho_inf,GAMMA);
/* q2_inf = 0.5*rho*U^2 */
q2_inf = 0.5*rho_inf*(GAMMA*pre_inf/rho_inf)*mach_inf*mach_inf;

fprintf(fpt,"Free stream values\n");
fprintf(fpt,"Mach              = %f\n",mach_inf);
fprintf(fpt,"Pressure          = %f\n",pre_inf);
fprintf(fpt,"Total Pressure    = %f\n",totpres_inf);
fprintf(fpt,"Temperature       = %f\n",temp_inf);
fprintf(fpt,"Total Temperature = %f\n",totemp_inf);
fprintf(fpt,"Entropy           = %f\n",entropy_inf);
fprintf(fpt,"\n");

for(i=0;i<no_spl_points;i++){
	point = spl_point[i];
	sonic = sqrt(GAMMA*p[point].var[4]/p[point].var[0]);
	rho = p[point].var[0];
	u = p[point].var[1];
	v = p[point].var[2];
	w = p[point].var[3];
	vel = sqrt(u*u + v*v + w*w);
	temp =  p[point].var[4]/(GAS_CONST*p[point].var[0]);
	pres = p[point].var[4];
	mach = vel/sonic;
	fact = 1.0 + 0.5*(GAMMA-1.0)*mach*mach;
	totemp = temp*fact;
	fact   = pow(fact,GAMMA/(GAMMA-1.0));
	totpres = pres*fact;
	entropy = pres/pow(rho,GAMMA);

	fprintf(fpt,"Point no. = %d\n",point);
	fprintf(fpt,"x = %f, y = %f, z = %f\n",p[point].x,p[point].y,p[point].z);
	fprintf(fpt,"Mach              = %f\n",mach);
	fprintf(fpt,"Pressure          = %f\n",pres);
	fprintf(fpt,"Total Pressure    = %f\n",totpres);
	fprintf(fpt,"Temperature       = %f\n",temp);
	fprintf(fpt,"Total Temperature = %f\n",totemp);
	fprintf(fpt,"Entropy           = %f\n",entropy);

	}
fclose(fpt);

return;
fpt = fopen("VG.OUT", "w");
fcp = fopen("CP.DAT", "w");
fprintf(fpt, "scalars  pressure\n");
for(i=0;i<no_grid_points;i++)
	if(p[i].y == 0.05){
		fprintf(fpt, "%18.6e\n", p[i].var[4]);
		cp = -(p[i].var[4] - pre_inf)/q2_inf;
		if(p[i].point_type==SOLID_WALL || p[i].point_type==KUTTA)
			fprintf(fcp, "%18.6e %18.6e\n", p[i].x, cp);
		}
fclose(fpt);
fclose(fcp);
}
/* -------------------------------------------------------------------------- */
void totalTemp(struct node* p){
/* -------------------------------------------------------------------------- */
int i;
double totemp,tomax,tomin,totemp_inf;
double sonic,vel2,mach2;
FILE* fp;

printf("Total temperature error into point.dat\n");

totemp_inf = temp_inf*( 1.0 + 0.5*(GAMMA-1.0)*mach_inf*mach_inf);

tomax = 0.0;
tomin = 10000000.0;

for(i=0;i<no_grid_points;i++){
    vel2 = p[i].var[1]*p[i].var[1] + p[i].var[2]*p[i].var[2] + p[i].var[3]*p[i].var[3];
    sonic = GAMMA*p[i].var[4]/p[i].var[0];

    mach2 = vel2/sonic;
    totemp = p[i].var[4]*(1.0 + 0.5*(GAMMA-1.0)*mach2)/(GAS_CONST*p[i].var[0]);
    tomax = max(tomax,totemp);
    tomin = min(tomin,totemp);
    }

fp = fopen(PointFile,"a");
fprintf(fp,"Percentage difference in total temperature = %f\n",(tomax-tomin)*100.0/totemp_inf);
fclose(fp);
}
/* -------------------------------------------------------------------------- */
double ERRF(double X){
/* -------------------------------------------------------------------------- */
double ARG,E,VB,T,ERR;

ARG = X*X;
E = 0.0;
if(ARG < 20.0) E = exp(-ARG);
VB = fabs(X);
T = 1./(1+0.3275911*VB);
ERR = 1.061405429*T;
ERR = (ERR-1.453152027)*T;
ERR = (ERR+1.421413741)*T;
ERR = (ERR-0.284496736)*T;
ERR = (ERR+0.254829592)*T;
ERR = 1-ERR*E;
if(X < 0.0) ERR = -ERR;
return ERR;
}
